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Abstract 

The  present  paper  concentrates  on  the  basic  equations  of  three-dimensional  problems  for 
nonlinear  piezoelectric  materials  of  hexagonal  systems  with  symmetry  class  6mm  and  of  cubic 
systems  with  symmetry  class  m3.  Emphasis  is  placed  on  developing  the  nonlinear  constitutive 
relations  between  extended  traction  (including  elastic  stress  and  polarization)  and  extended 
strain  (including  elastic  strain  and  electric  held).  The  corresponding  one-dimensional 
mathematical  models  for  piezoelectric  ceramic  with  symmetry  classes  6 mm  and  m3  are  also 
given.  Numerical  examples  are  also  carried  out  for  the  impact  problem  to  show  the  important 
effect  of  the  piezoelectric  nonlinearity  on  the  stress  wave.  Therefore,  the  derived  concise 
equations  can  directly  be  applied  to  evaluate  the  nonlinear  piezoelectric  effects  of 
piezoelectricity  by  the  nonlinear  finite  element  method. 


1.  Introduction 

The  phenomenon  of  piezoelectricity  was  first  discovered 
by  the  Curie  brothers,  Pierre  and  Jacques,  in  1880. 
Piezoelectricity  is  currently  enjoying  a  great  resurgence  in  both 
fundamental  research  and  technical  applications  (Chen  1971, 
Chizhikov  et  al  1982,  Maugin  1988,  Ashida  and  Tauchert 
1998,  Chandrasekharaiah  1998).  Although  the  behavior  of 
piezoelectric  materials  in  non-structural  applications  has  been 
investigated  extensively,  the  treatment  is  often  simplistic.  In 
particular,  most  of  the  achievements  have  only  been  made  in 
the  frame  of  linear  constitutive  relations. 

The  rapid  development  of  computer  science  and  nonlinear 
finite  element  applications  reveals  the  importance  of  estab¬ 
lishing  the  nonlinear  basic  theory  for  piezoelectricity.  Nelson 
(1978),  Toupin  (1983),  Tiersten  (1981)  studied  the  nonlinear 
theory  of  dielectrics.  Norwood  et  al  (1991),  Kulkami  and 
Hanagud  (1991)  used  a  Neo-Hookean  constitutive  relation 
to  model  the  response  of  piezoelectric  ceramics.  Pai  et  al 
(1992)  considered  the  dependence  of  the  piezoelectric  strain 
parameters  upon  the  strain  in  formulating  a  plate  theory  of 
piezoelectric  laminates.  Joshi  (1992)  considered  the  nonlinear 
constitutive  relations  for  piezoelectric  materials,  where  a 
concise  expression  was  given.  Tiersten  (1993)  investigated 
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the  nonlinear  problems  of  thin  plates  subjected  to  large  driving 
voltages.  However,  most  of  these  studies  considered  only  the 
case  of  small  deformations,  and  second-order  items  were  often 
neglected.  Later  on,  based  on  the  theory  of  invariants,  from 
invariant  polynomial  constitutive  relations,  Yang  and  Batra 
(1995)  investigated  the  second-order  theory  for  piezoelectric 
materials  with  symmetry  class  6 mm  and  class  mm2,  where 
only  nine  independent  stiffness  constants  were  introduced  for 
the  nonlinear  items. 

In  this  paper,  we  develop  the  basic  constitutive  equations 
for  three-dimensional  nonlinear  problems  of  piezoelectric 
materials.  The  polynomial  constitutive  relations  for  the 
6 mm  crystal  class  and  m3  crystal  class  are  derived  using 
an  invariant  integrity  basis.  Furthermore  all  the  basic 
equations  are  specialized  to  those  corresponding  to  the  one¬ 
dimensional  model.  This  one-dimensional  model  is  then 
solved  by  implementing  the  governing  equations  in  the 
COMSOL  Multiphysics  software  (COMSOL  2008).  These 
numerical  examples  of  the  impact  problem  show  the  influence 
of  the  piezoelectric  nonlinearity  on  the  stress  wave.  Since 
the  present  theory  is  not  restricted  to  small  applied  loads 
or  electric  fields,  it  should  therefore  be  very  useful  to 
researchers  for  the  investigation  of  the  mechanics  and  physics 
of  piezoelectricity  undergoing  large  nonlinear  deformations 
(Cheng  1996,  Hokstad  2004). 
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2.  Statement  of  dynamic  problems  for  a  nonlinear 
piezoelectric  material 


2.2.  Boundary  and  initial  conditions 

(a)  Mechanical  boundary  conditions 


Let  the  coordinates  of  a  material  particle  with  respect  to 
a  rectangular  Cartesian  coordinate  system  be  Xk  in  the 
reference  configuration  (undeformed  configuration)  and  its 
spatial  coordinates  in  the  current  configuration  (deformed 
configuration)  be  x*. 


2.7.  Governing  equations 

For  a  nonlinear  piezoelectric  material,  in  the  absence  of 
body  force  and  free  charge  density,  the  moving  equation  and 
the  quasistatic  approximation  to  Maxwell’s  equations  in  the 
reference  configuration  can  be  written  as  (Yang  and  Batra 
1995,  Pao  1978,  Kiral  and  Eringen  1990) 

[ TklX/c.l  +  J  F 'Kl  Sq(F ’m\F 'nI  —  \FMltFN^8u)EMEN]tK 


=  PoSlkUk, 

(1) 

TIk.k  +  J  soF Kl(F m\E  m)  ,k  —  0, 
where 

(2) 

J  =  det  (F) , 

9x 

(3) 

IS 

II 

(4) 

Tkl  is  the  second  Piola-Kirchhoff  (PK2)  stress  tensor,  u  K 
the  mechanical  displacement  vector,  n  K  the  material  electric 
polarization,  Ek  the  electric  field,  po  the  mass  density  in 
the  reference  configuration,  £o  the  permittivity  of  free  space, 
SicK  the  shifter,  and  S u  the  Kronecker  delta.  Throughout  this 
paper,  a  repeated  index  implies  summation  over  the  range 
of  the  index,  and  a  comma  followed  by  K(i)  implies  partial 
differentiation  with  respect  to  Xk(x/).  A  dot  above  a  quantity 
signifies  its  material  time  derivative.  In  addition,  the  PK2  stress 
tensor  Tkl  can  be  written  as 


Tkl  =  JXK,kXL,i  {aa  +  PkE,) ,  (5) 

where  I),  and  El  are,  respectively,  the  Cauchy 

stress,  electric  polarization,  and  electric  field  in  the  current 
configuration. 

For  an  isothermal  or  adiabatic  process,  in  the  material 
configuration,  the  energy  density  function  in  the  absence  of 
heat  conduction  or  heat  sources  can  be  defined  as  (Pao  1978) 

X  =  X  (T KL,  Ek)  ,  (6) 


where  Tkl  is  the  Green-Lagrange  strain  tensor,  and  E  K  the 
electric  field  in  the  material  form. 

Also,  the  displacement-strain  and  electric  potential- 
electric  field  can  be  expressed  as 

T  KL  =  \  (xk,KXk,L  —  &Kl)  =  \  ( liK.L  +  UL.K  +  UM.RUM.l)  , 

(7) 

Ek  =  —<P,k-  (8) 

In  addition,  the  basic  equations  (1)  and  (2)  are 
accompanied  by  the  following  constitutive  relations  (Pao  1978) 


Tkl  = 


nK  = 


9X 

drKL’ 

9E 

~dE^ 


(9) 

(10) 


uk  —  uk  X  g  Si,  (11) 

u°kxi,l  (Tkl  —  PoJ  1TIkEl)  Sim  =  Pm  X  e  S2, 

(12) 

where  5)  U  S2  =  S,  Si  fl  S2  =  0,  and  S  covers  the 
total  boundary;  uk  and  pM  are,  respectively,  the  given 
displacement  and  traction  on  the  boundary  in  the  reference 
configuration;  n°K(X)  is  the  unit  normal  to  the  body  in  the 
reference  configuration. 

(b)  Electrical  boundary  conditions 

0  =  0  Xe  S3,  (13) 

n°K  (7  x PqTIk  +  £qEk)  —  Dno  X  e  54,  (14) 

where  S3  U  S4  =  S,  S3  fl  S4  =  0;  0  and  Dno  are, 
respectively,  the  known  electric  potential  and  normal 
electric  displacement  on  the  boundary  in  the  reference 
configuration. 

(c)  Initial  conditions 

uk(X,0)  =  u°k(X),  (15) 

uK  (X,  0)  =  u°K  (X)  ,  (16) 

where  u°K  is  the  known  initial  displacement  and  u°K  the 
known  initial  velocity. 

3.  Transversely  isotropic  materials  with  symmetry 
class  6 mm 

3.1.  Polynomial  integrity  basis 

Two  basic  requirements  of  invariance  that  must  be  imposed 
upon  the  constitutive  equations  are  spatial  invariance  and 
material  invariance  (Jordan  and  Eringen  1964).  For  6 mm 
piezoelectricity  with  both  the  elastically  symmetric  axis  and 
the  poling  direction  as  the  SVaxis,  the  polynomial  integrity 
basis,  the  degree  of  which  is  less  than  three,  can  be  given  in  the 
following  concise  form  (Kiral  and  Eringen  1990). 

Elements  in  T/y  only: 


Degree  1: 

r33,  Hi  +  r22, 

(17) 

Degree  2: 

rnr22  -  t\2,  ri3  +  r23> 

(18) 

Degree  3: 

r„  (r?!  +  6r„r22  -  i2r?2  +  9Fj2)  , 
r  1 1 1^23  +  r22r^3  —  ri3r23ri2. 

(19) 

Elements  in  E;  only: 

Degree  1:  E  3, 

(20) 

Degree  2:  E\  +  E\. 

(21) 

Elements  in  F/y  and  E /  only: 

Degree  2:  r3i£i  +  T23E2, 

(22) 

Degree  3: 

(£iF23  +  £2r31)ri2-£1r22r31-£2r 

1 1  r23 , 

(23a) 

Fh£2  +  T21E-l  —  2EiE2T  i2. 

(23  b) 
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3.2.  Polynomial  free  energy  function 

In  order  to  derive  the  second-order  nonlinear  constitutive 
equations,  the  energy  density  function  E  should  be  formed  as 
a  third-order  polynomial  function  of  T/y  and  E / .  Neglecting 
the  initial  stress  and  initial  polarization,  i.e.  assuming  that  the 
energy  density  function  of  the  piezoelectricity  is  zero  when  it 
is  in  a  free  state,  from  equations  (17)  to  (23),  after  complicated 
symbolic  mathematical  manipulations  using  the  Mathematica 
software  package,  the  energy  density  function  E  can  finally  be 


obtained  in  the  following  form 

S  =  E!„  +  E£  +  Ee  +  Ec  +  +  Eg  +  Eg,  (24) 

where 

=  a\  (rur22  —  r22)  +  «2  (rf3  +  r31)  +  a3r23 

+  a4  (r  1 1  +  r22)  r33  +  <25  (F 1 1  +  r22)2  ,  (25a) 

E£  =  e,  (Ej  +  Ej)  +  e2E2,  (25 b) 

Ef  =  e\  (r3i£i  +  T23E2)  +  e2r33 E3 

+  e3  (r  1 1  +  r22)  £3 .  (25c) 


Ec  =  c,r n  (r2;  -  I2r22  +  6r„r22  +  9r22) 

+  c2  (r  1 1  rf3  +  r3i  (r22r3i  —  ri2r23)) 

+  c3  (r  11  r22  —  r22)  r33  +  C4  (rj3  +  r31 )  r33 
+  C5  (rn  +  r22)  (rnr22  —  r22) 

+  Ce  (r  11  +  r22)  (r23  +  r31)  +  c7r33 

+  c8  (r  1 1  +  r22)  r33  +  Cg  (rn  +  r22)3 

+  C10  (rn  +  r22)2  r33,  (25  d) 

E,  =  rn  [E\  +  El)  E3  +  mE\,  (25e) 

Sg  =  g\  (ri2r23.Ei  —  r22r31£1  —  rur23,E2  +  r \2t3\E2) 

+  g2r33  (r3i£i  +  r23£2) 

+  g3  (Tn  +  r22)  (r31£j  +  r23£,2) 

+  ^4  (— rj2  +  rur22)  e3  +  g5  (r23  +  r31)  e3 
+  ^r33£3  +  gi  (r  1 1  +  r22)2  e3 
+  gs  (r  1 1  +  r22)  r33£3,  (25/) 

SG  =  Q\  (r 21e\  -  2F12£1£2  +  TnE\) 

+  Q2r33  (e\  +  Ej)  +  Q3  (r  1 1  +  r22) 
x  (E\  +  Ej)  +  Q4  (r3X£i  +  r23E2)  e3 
+  QsE33eI  +  q6  (rn  +  r22)  e\.  (25 g) 

Equation  (24)  indicates  that  the  energy  density  consists 
of  seven  parts.  Obviously,  E a,  Ee,  and  E(,  correspond  to 
the  linear  constitutive  relations,  and  Ec,  E,,  Eg,  and  E q 
correspond  to  the  nonlinear  constitutive  relations.  Also,  it  is 
easily  observed  that  there  are  a  total  of  36  independent  material 
constants  for  the  6 mm  crystal  class,  which  is  in  agreement 
with  the  early  literature  (Landolt  1966),  and  that  there  are  ten 
independent  constants  for  the  corresponding  linear  material, 
which  is  further  in  agreement  with  the  well-known  results. 

3.3.  Second-order  constitutive  theory 

Substituting  equation  (24)  together  with  equation  (25)  into 
equations  (9)  and  (10),  we  obtain  the  elastic  stress  and  electric 


polarization  as  follows 

Tij  =  TUa  +  TIJe  +  Tjjc  +  Tjjg  +  Tijq,  (26) 

11/  =  Il/E  +  11  [e  +  I!/,;  +  Il/g  +  TIiq,  (27) 

where 

T\\a  =  air22  +  a4r33  +  2a$  (r u  +  r22) ,  (28a) 

T\\e  =  e3E3,  (28  b) 

Tuc  =  3 Ci  (r2!  -  4r22  +  4r„r22  +  3r22) 

+  C2r23  +  c3r22r33  +  c3  (2r nr22  +  r22  —  r^) 

+  c6  (r23  +  r2))  +  c8r23  -  3  c9  (r„  +  r22)2 


+ 

2C10  (r  11  +  r22)  r33, 

(28  c) 

T\\g 

—gi  r23 e2  +  g3  (r3i£'i  +  t23e2) 

■ 

g4E22E3  +  2g7  (ru  +  r22)  e3  +  gsr33E3, 

(28  d) 

T\\q 

-QxE22  +  Q3  (El  +  El)  +  Q6E23, 

(28c) 

T22a 

aiTn  +  a4r33  +  2 a3  (Fn  +  r22) , 

(29a) 

T22e 

e3E3, 

(29  b) 

T22C 

-eCiTn  (r„  +  3r22)  +  c2r2, 

■ 

c3r„r33  +  c5  (2r„r22  +  r2!  -  r22) 

+ 

c6  (r23  +  r2))  +  c8r23  +  3 c9  (rn  +  r22)2 

+ 

2C10  (ru  +  r22)  r33, 

(29  c) 

T22g 

— ^ir3i£i  +  g3  (r3i£'i  +  r23£2) 

• 

^rn£3  +  2g7  (ru  +  r22)  e3  +  g$r33E3, 

(29  d) 

T22q 

-  Q\E\  +  Q3  (El  +  E2)+Q6E2, 

(29c) 

T33a 

2a3F33  +  04  (ru  +  r22) , 

(30a) 

T33e 

e2E3, 

(30fe) 

T33c 

■■  c3  (r„r22  -  r22)  +  c4  (r23  +  r32,) 

■ 

3C7r23  +  2C8  (r„  +  r22)  r33 

+ 

C10  (r  11  +  r22)2 , 

(30c) 

T33g 

= 

g2  (r3i  £■  1  +  r23  e2)  +  2g6r33E3 

+ 

g&  (Fn  +  r22)  e3, 

(30rf) 

T33q 

Q2  (El  +  El)  +  Q5E23, 

(30c) 

T23a 

a2F23, 

(31a) 

T23e 

0.5ei£2, 

(31/7) 

T23c 

;  c2(r„r23  —  0.5  r  i2r3i) 

■ 

c4r23r33  +  C6  (Fn  +  r22)  r23, 

(31c) 

T23g 

o.5gi  (r \2E\  —  r  u  e2)  +  o.5g2r33.E2 

• 

o.5g3  (ru  +  r22)  e2  +  g5r23E3, 

(31  d) 

T23d 

= 

0.5  q4e2e3. 

(31c) 

T3\a 

a2r3i, 

(32a) 

T3\e 

O.SeiTii, 

(32/7) 

T3ic 

=  c2  (r22r3i  —  o.5r  12r23)  +  c4r31r33 

• 

Q  (Fu  +  r22)  r3i, 

(32c) 

T3ig 

o.5gi  (r  1  2e2  —  r22Ei)  +  o.5g2r33£, 

■ 

o.5g3  (ru  +  r22)  E\  +  g5r3i£,3, 

(32d) 

T3\q 

1  Q.5Q4E\E3, 

(32c) 
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Tua  = —a\V\2,  (33a) 

Tnc  =  —  i2CiFnri2  —  o.5C2r23r3i 

—  c3r12r33  -  c4r12  (ru  +  r12) ,  (33  b) 

Tne  =  0,  (33c) 

Tng  =  0.5^1  (r23£j  +  r3i£2) ,  (33 d) 

Ti2Q  =  —  Q\E\E2,  (33e) 

IllE  =  —2e\Ei,  (34  a) 

Ylle  =  -eir31,  (34ft) 

nh;  = -2?7  \EiE3,  (34  c) 

riig  =  —  gi  (r i2r23  —  r22r3i)  —  g2r3iF33 

—  gs  (r  1 1  +  r22)  r3i,  (34c/) 

riig  =  2Q1  (r22£j  —  r  32£2)  —  2Q2r33£,i 

—  203  (T j i  +  P22)  E\  —  Q4r3\E3,  (34e) 

n2£  =  -2  elE2,  (35  a) 

n2e  =  —  eir23,  (35  b) 

n2ri  = -2rnE2E3,  (35c) 

n2g  =  —  gi  (r  i2r3i  —  r  nr23)  —  g2r23r33 

—  g3  (Fn  +  r22)  r23,  (35 d) 

n2G  =  2Qi  (r 1 1 £2  —  r i 2£i)  —  2Q2r33£,2 

—  2(?3  (T 1 1  +  r22)  E2  —  Q4r23E3,  (35e) 

n3s  =  -2e2E3,  (36  a) 

n3e  =  — e2r33  —  e3  (r n  +  r22) ,  (36 b) 

n3,  =  -rn  (Ef  +  El)  -  3/72 El  (36c) 

n3j?  =  —gi  (rnr22  -  r22)  -  g5  (r23  +  r3l ) 

—  —  gi  (Tn  +  r22y  —  g%  (ru  +  r22)  r33,  (36 d) 
n3G  =  —Qa  (r 31  -Ei  +  r23£2)  —  2Q3y33e3 

-2Q(s{Tn  +  r22)  E3.  (36e) 


Equations  (28)-(36)  can  easily  be  written  in  matrix  forms. 
Thus,  we  have  obtained  the  concise  expressions  of  nonlinear 
constitutive  equations  for  the  6mm  crystal  class,  which  is  the 
basis  for  analyzing  the  nonlinear  problems  of  piezoelectric 
materials  with  this  kind  of  symmetry. 

It  should  also  be  pointed  out  that  the  corresponding  results 
of  small  deformations  and  weak  electric  fields  (Yang  and  Batra 
1995)  and  the  corresponding  results  of  small  deformations  and 
strong  fields  (Tiersten  1993,  Yang  and  Batra  1995)  can  easily 
be  deduced  from  equations  (26)-(36). 

4.  Cubic  materials  with  symmetry  class  m3 

For  m3  piezoelectricity,  the  polynomial  integrity  basis,  the 
degree  of  which  is  less  than  three,  is  given  by  (Kiral  and 
Eringen  1990). 

Elements  in  F/y  only: 

(37) 


Degree  3:  rnr22r33,  r23r3iri2, 

(39a) 

r  IlTjj  +  F  2^3!  +  F33rj'2, 

rn  (r3i  -  rr2)  +  r22  (r\2  -  f223)  +  r33  (r223  -  r2,) ,  (39 b) 
Ti  i  r22  (rj  i  —  r22)  +  r22r33  (r22  —  r33) 

+  r33Fn  (r33  —  Tn) .  (39c) 

Elements  in  Ej  only: 

Degree  2:  E2{  +  £2  +  is2.  (40) 

Elements  in  F/y  and  E /  only: 

Degree  3:  £’2£’3r23  +  E3E\Y3\  +  E1E2Y i2, 

£2r„  +  £2r22  +  £2r33,  (41a) 

El  (r22  -  r33)  +  e\  (F33  -  r„)  +  e\  (r„  -  r22) .  (4\b) 

Similarly,  the  polynomial  energy  density  function  can  be 
expressed  as 

X  =  Xfl  +  X£  +  XG  +  XG,  (42) 

where 

xn  =  a\  (rn  +  r22  +  r33)2 

+  a2  (r22r33  +  r33rn  +  rnr22) 

+  a3  (F23  +  r2j  +  r22)  ,  (43a) 

Sc  =  Cl  (Fn  +  r22  +  r33)3  +  c2  (ru  +  r22  +  r33) 
x  (r22r33  +  r33rn  +  rnr22) 

+  c3  (rn  +  r22  +  r33)  (r23  +  r31  +  r22) 

+  C4rnr22r33  +  Csr  i2r23r3i 
+  C(,  (r 1 1 rf3  +  r22r31  +  r33ri'2) 

+  c7(rf2r22  —  r23r22  +  r31r  n  —  rj"2r  3 1  +  r23r33 
—  r31r33)  +  (r n  —  r22)  (r n  —  r33)  (r22  —  r33) , 

(43b) 

X£  =  £l  (El  +  Ej  +  El) ,  (43c) 

SG  =  Qx  (r„  +  r22  +  r33)  (e\  +  e\  +  e]) 

+  G2  (E12E1E2  +  r3iExE3  +  F 23E2E2,) 

+  G3((F22  -  r33)  El  +  (F33  -  r n)E\ 

+  (r„  -  r22)£2)  +  q4  (rnEl  +  y22e22  +  r33£2) . 

(43  d) 

Substituting  equations  (42)  and  (43)  into  equations  (9) 
and  (10),  we  finally  obtain  the  elastic  stress  and  electric 
polarization  as  follows 

Tjj  =  Tija  +  Tuc  +  Tijq,  (44) 


Degree  1: 
Degree  2: 


Tn  +  r22  +  r33, 
rnr22  +  r22r33  +  r33rn, 


(38) 


ii/  =  ii/E  +  ii/G, 


(45) 
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c0H(t) 


c0H(t) 


Figure  1.  A  one-dimensional  piezoelectric  model  of  6 mm  crystal 
class  loaded  at  the  left  end  along  the  poling  direction. 


Figure  2.  A  one-dimensional  piezoelectric  model  of  6 mm  crystal 
class  loaded  at  the  left  end  normal  to  the  poling  direction. 


where 

T\\a  =  2«i  (Tn  +  r22  +  r 33)  +  a2  (r22  +  r33)  ,  (46a) 

t’i  ic  =  3Ci  (r  ii  +  r22  +  r  33)“ 

+  c2((r22  +  r33)  (Tn  +  r22  +  P33)  +  rnr22 
+  r  1 1  r33  +  r22r33)  +  c3  (rf2  +  rf3  +  r31) 

+  c4r22r33  +  +  Ci  (rjj  —  r^) 


+  c8  (2r„r22  -  2r„r33  +  r323  -  r22) ,  (46*) 

Tuq  =  01  (Ef  +  E\  +  £2)  +  03  (£3  -  El)  +  Q4El 

(46  c) 

r23a  =  03^3,  (47a) 

Tuc  =  c3r23  (r  1 1  +  r22  +  r33)  +  o.5C5F3iri2 

+  C6rnr23  +  c2r23  (r33  —  r22) ,  (47*) 

T23Q  =  O-502£2£3,  (47c) 

FIie  =  —  2ei£i,  (48a) 

n  1  q  =  —20i  (r„  +  r22  +  r33)  £1  —  02  (r  i2£2  +  r3i£3) 
-  203  (r22  -  r33)  £1  -  204ru£i.  (48*) 


The  other  stress  and  polarization  components  can  be 
obtained  by  simple  index  permutations,  which  are  omitted  here 
for  brevity. 

As  shown  in  equations  (46)-(48),  for  linear  cubic 
piezoelectricity,  there  are  five  independent  material  constants, 
which  is  in  agreement  with  the  well-known  results.  It  is 
interesting  to  note  that  for  the  nonlinear  m3  crystal  class,  the 
squared  terms  of  strain  and/or  electric  field  components  have 
no  effects  on  the  electric  polarization,  whilst  the  products  of 
strain  and  electric  field  components  have  no  contribution  to 
the  stress  distribution.  In  addition,  the  analysis  on  the  cubic 
system  with  m3  crystal  class  also  shows  that  even  for  nonlinear 
materials,  there  are  only  a  total  of  16  independent  material 
constants. 

Furthermore,  we  remark  that  if  C7,  C8,  and  03  are  all 
set  to  be  zero,  then  equations  (46)-(48)  directly  reduce  to  the 
corresponding  constitutive  equations  of  the  cubic  system  with 
the  symmetry  class  m3m. 


5.1.  Mechanical  loads  parallel  to  the  poling  direction 

In  this  case,  as  shown  in  figure  1,  =  X3  +  w(X 3,  t )  = 

Z  +  w(Z,  0,  M3  =  w(Z,t),  4>  =  4>(Z,t),  the  constitutive 
equations  can  be  directly  simplified  from  the  three-dimensional 
ones  as 


T33  =  ^3303  —  £33 £3  +  ^0333 £33  —  ^  0333 £3  ~  <?333r33£3> 

(49) 

n3  =  e33T33  +  £33  £3  +  ^g333r23  +  ^333  £3  +  0333£33£3, 

(50) 

where  a3 3,  633,  £33,  C333,  £333,  0333,  and  77333  are  the  new 
and  independent  material  constants,  which  are  introduced  for 
clarity.  Also, 


r33  (Z,  0  = 


3  w  1  /  3  w 
3 Z  +  2I 3Z 


£3  (Z,  t)  = 


dtp 
3 Z' 


(51) 

(52) 


The  governing  equations  can  finally  be  expressed  as 


duA 


d2w 


T33, z  I  1  +  —  I  +  Th——  +  £0£3£3,z  1  +  —  ) 


£o  £3 


3Z  J 
,d2w 
3Z2 


1  + 


3Z2 
Aw 
3 Z 


dw\ " 


3  Z 


-3 


p0w, 


(53) 


n3,z  —  so  1 1  + 


Aw 

dZ 


1  ,  dw 

S»l1  +  dZ 


2  d2w 

Zz2' 


£3,z  =  0. 


(54) 


For  the  mechanical  impact  only,  the  mechanical  and 
electric  boundary  conditions  are  given  as 

du)  (0,  t )' 


T33  (0,  t)  ^1  + 

=  -o0H  (t) , 
w  (l,  t )  =  0, 


dZ 


-  po  (0)  n3  (o,  t)  £3  (0,  t) 


(55  a) 
(55*) 


5.  One-dimensional  models  of  nonlinear  impacts  for 
the  6 mm  crystal  class 

As  shown  in  figures  1  and  2,  a  bar  of  rectangular  cross-section 
with  length  l  is  fixed  at  the  right  end.  In  this  section,  we 
will  derive  the  simplified  equations  of  one-dimensional  impact 
problems  under  two  kinds  of  mechanical  impact  loadings. 


Aw  (0,  t)\ 

1  +  dz  j  Po  (0)  10  (0,  t)  +  £0£3  (0,  t )  =  0,  (56a) 

Aw  (l,  t )  \ 

1  +  dz  j  Po  (/)  n3  (1,0  +  £0£3  (/,  0  =  0,  (56*) 
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where  ap  is  the  known  constant  and  H(t)  the  Heaviside  unit 
step  function. 
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Table  1.  Nonlinear  piezoelectric  coefficients  for  quartz. 


Name 

Present 

D&Ga 

Value 

Second-order  elastic  stiffness  (GPa) 

a33 

CE 

'-33 

8.6736  x  10lob 

Piezoelectric  coupling  (C  m~2) 

633 

£33 

0. 17 lc 

Dielectric  constant  (eod) 

£33 

b33 

4.40c 

Third-order  elastic  (GPa) 

C333 

CE 

*-333 

-30.0  x  10lof 

Third-order  piezoelectric  (£°  m  V-1) 

#333 

1  jklm 

—  1.31s 

Electrostrictive  (e°  F  m~') 

<2333 

/333 

— 4.40h 

Third-order  dielectric  (F  V-1) 

7)333 

p "1 

fc333 

0(— 3.5  x  10-17)1 

11  D  &  G*:  Davison  and  Graham  (1979).  b  Value  for  quartz  from  COMSOL  3.5a 
Material  Library.  c  Value  from  D&G. 
d  e°  =  8.854  x  1(L12  F  m_1  (electric  constant  in  vacuum). 
e  Relative  dielectric  constant  value,  references  to  Bechmann  (1958)  and  Fontanella 
and  Andeen  (1974). 

f  Same  order  as  the  (second-order)  elastic  stiffness  but  significantly  large  for  finite 
strain,  reference  to  D&G. 

g  The  contribution  is  typically  several  orders  smaller  than  other  terms,  reference  to 
D&G. 

h  The  electrostrictive  coefficient  is  typically  the  same  order  as  the  second-order 
dielectric  but  of  opposite  sign. 

1  Several  orders  smaller  than  other  coefficients  and  can  be  neglected  unless  the 
.E-field  is  extremely  high,  e.g.  several  orders  larger  than  the  strain. 


The  initial  conditions  can  be  expressed  in  terms  of  the 
known  displacement  w°(Z)  and  velocity  w°(Z)  as  follows 

w(Z,  0)  =  w°(Z),  (57a) 

ii>(Z,  0)  =  w°(Z).  ( 57b ) 

It  should  be  noted  that,  as  shown  in  equations  (49) 
and  (50),  when  the  loading  direction  is  in  the  poling  direction, 
there  are  a  total  of  seven  independent  material  constants,  which 
include  three  linear  coefficients  (i.e.  a 33,  633,  and  £33)  and  four 
nonlinear  coefficients  (i.e.  C333,  £333,  Q333,  and  77333). 


5.2.  Mechanical  loads  vertical  to  the  poling  direction 


In  this  case,  as  shown  in  figure  2,  x2  =  Y  +  v(Y,  t ),  M2  = 
v(Y,  t),  cp  —  <p(Y,  t ),  the  constitutive  equations  are  obtained  as 


T22  =  a  iiT22  +  \C222T222  —  \Q222E2,  (58) 

n2  =  SUE2  +  Q222Y22E2,  (59) 


where  an,  e n,  C222,  and  Q 222  are  the  new  material  constants 
introduced  here  for  clarity.  Also, 


r22  (Y,  0  = 


dv 

8Y 


(60) 


Correspondingly,  the  conditions  of  unique  solution  can  be 
described  as 

/  du  (0,  t)\ 

t22  (0,  t)  h  +  ^y  j  -  po  (0)  n2  (0,  t )  e2  (0,  t ) 


=  —(JqH  (t) ,  (64a) 

v  (l,  t )  =  0,  (64b) 

l  du  (0,  t)\ 

(  1  +  ^y  J  Po  (0)  n2  (0,  t )  +  e0E2  (0,  t )  =  0,  (65a) 

(  d  v(l,t)\ 

(  1  +  — j  Po  (0  n2  (Z,  0  +  s0E2  (l,t)  =  0,  (65b) 

v  (Y,  0)  =  u°  (T) ,  (66a) 

v  (Y,  0)  =  v°  (Y) ,  (66 b) 


where  v°(Y)  and  v°(Y)  are,  respectively,  the  known  initial 
displacement  and  velocity. 

As  shown  in  equations  (58)  and  (59),  when  the  loading 
direction  is  vertical  to  the  poling  direction,  there  are  a  total  of 
four  independent  material  constants,  which  include  two  linear 
coefficients  (i.e.  an  and  £n)  and  two  nonlinear  coefficients 
(i.e.  C222  and  Q222). 


E2  (Y,  t)  =  - 


dY 


T22l] 


du  \ 

1+a?) 


T22-^^  +  eaE2E2j 


du 

1  +  9? 


d2v 

—  £0  Ej r 

2  dY2 


1  + 


di> 

8Y 


-3 


Pqv, 


n 


2  ,Y 


SO  (  1  + 


+  SO 


du 

dY 


du 
1  +  dY 


d2v 

d?2 


E  2 


E2,y  =  0. 


(61) 


The  governing  equations  can  finally  be  expressed  as 

d2v  2  J-x_2 


(62) 


(63) 


6.  One-dimensional  model  of  nonlinear  impacts  for 
the  m3  crystal  class 

For  the  one-dimensional  m3  crystal  class  structure,  if  the 
applied  mechanical  loads  are  parallel  to  the  X2(Z)  axis  (see 
figure  3),  we  have  x2  =  2^3  +  w(X2,t)  —  Z  +  w(Z,t), 
m3  =  w(Z,  t),  (p  =  (p(Z,  t). 

Similar  to  the  6 mm  crystal  class,  the  governing 
equations  (53)  and  (54),  the  extended  geometry  equations  (51) 
and  (52),  and  the  conditions  of  unique  solution  (55)— (57)  all 
remain  the  same.  The  constitutive  equations  can  be  expressed 
as 

Tl3  =  ^33^33  +  5C333F323  ~  \Q333E3,  (67) 
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Figure  3.  A  one-dimensional  piezoelectric  model  of  m3  crystal  class 
loaded  at  the  left  end. 


n3  =  S33E3  +  Q33sT33E3.  (68) 

As  shown  in  equations  (58)  and  (59),  for  the  one¬ 
dimensional  m  3  crystal  class  structure,  there  are  a  total  of 
four  independent  material  constants,  which  include  two  linear 
coefficients  (<233  and  £33)  and  two  nonlinear  coefficients  (C333 
and  g333)- 

7.  Numerical  results 

From  the  governing  equations,  equations  (1)  and  (2),  it  is  seen 
that  for  a  one-dimensional  response,  the  only  change  between 
crystal  classes  is  the  non-zero  constants  in  the  constitutive 
law  which  are  different  for  different  crystal  classes.  This  has 
been  explicitly  shown  for  the  6 mm  and  3m  crystal  classes 
in  section  6.  We  point  out  that  experimental  values  for  the 
nonlinear  coefficients  are,  in  general,  not  readily  available  for 
most  materials.  However,  values  for  quartz,  which  is  in  the 
32  crystal  class,  are  available  (Davison  and  Graham  1979). 
Consequently,  for  the  illustrative  purpose,  equations  (53) 
and  (54)  were  solved  numerically  using  quartz  properties. 
Although,  in  general,  the  coefficients  in  equations  (28)- 
(36),  even  for  small  deformations,  are  not  the  same  as  the 
classical  constitutive  coefficients  (Feng  et  al  2009),  in  the  one¬ 
dimensional  case,  however,  there  is  a  direct  correspondence. 
This  correspondence  between  those  adopted  in  this  paper  and 
those  in  Davison  and  Graham  (1979)  is  shown  in  table  1. 
These  values  were  used  when  solving  equations  (53)  and  (54) 
using  the  COMSOL  solver  (COMSOL  2008)  via  the  nonlinear 
constitutive  relations  (49)  and  (50).  To  reduce  numerical  noise, 
Rayleigh  damping  was  added  via  an  additional  weak  form 
term.  In  the  numerical  solutions,  three  different  models  are 
considered:  (1)  the  Dubner-Abate-Crump  (DAC)  model  for 
the  linear  small  strain  solved  via  Laplace  transform  techniques 
using  a  modified  DAC  algorithm  (Laverty  and  Gazonas  2005); 
(2)  the  finite  strain  model  where  the  finite  strain  terms  in  the 
governing  equations  and  constitutive  relations  equation  (49)— 
(5 1)  are  included  in  the  solution  but  nonlinear  effects  due  to  the 
electric  field  are  ignored;  (3)  the  full  nonlinear  model  where  the 
fully  nonlinear  governing  equations  (53)  and  (54),  along  with 
the  nonlinear  constitutive  relations  (equations  (49)  and  (50)), 
are  solved  using  the  corresponding  material  properties  in 
table  1 .  As  mentioned,  the  solutions  for  the  second  and  third 
models  were  obtained  by  implementing  the  relevant  equations 
in  COMSOL  (COMSOL  2008). 

First,  a  step  pressure  load  of  op  =  1  GPa  is  applied  at 
the  left  end  of  the  piezoelectric  quartz  bar  and  the  stress  wave 


Figure  4.  One-dimensional  stress  wave  response  in  the  midpoint  to  a 
step  pressure  loading  in  the  nonlinear  piezoelectric  quartz  bar:  a  load 
of  cr0  =  1  Pa  based  on  the  DAC,  linear  strain,  and  finite  strain  models 
in  (a),  and  a  load  of  op  =  10  GPa  based  on  the  DAC,  finite  strain,  and 
full  nonlinear  models  in  (b). 

(This  figure  is  in  colour  only  in  the  electronic  version) 


response  at  the  middle  point  (1/2)  is  calculated  (figure  1)  at 
different  times.  The  results  for  the  three  models  are  shown 
in  figure  4(a).  It  is  observed  that,  for  the  infinitesimal  strain 
deformation,  the  numerical  result  via  COMSOL  is  identical  to 
the  analytical  solution  using  the  modified  DAC.  This  indicates 
that  the  addition  of  the  Rayleigh  damping  does  not  adversely 
affect  the  solution.  It  is  also  noted  that  the  solution  based  on 
the  finite  strain  model  via  the  COMSOL  solver  is  also  very 
close  to  the  linear  strain  modified  DAC  solution  although  there 
is  a  slight  phase  shift  in  the  response.  This  clearly  illustrates 
that  quartz  continues  to  respond  linearly  up  to  moderately  large 
pressures. 

Second,  the  pressure  load  is  increased  to  op  =  10  GPa, 
which  corresponds  to  a  pressure  value  that  is  typical  of  shock 
loading,  and  the  stress  wave  responses  in  the  middle  point 
of  the  bar  are  calculated  (figure  4(b))  for  the  three  models. 
It  is  seen  that  at  this  loading  level  the  response  exhibits  a 
significant  phase  shift,  due  to  the  coupling  and  nonlinearity  in 
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the  constitutive  relations  (Lysne  1972).  It  is  further  observed 
that  at  a  high  loading,  the  amplitude  of  the  stress  wave  can  also 
be  significantly  increased  due  to  the  piezoelectric  nonlinearity. 
It  is  further  interesting  to  observe  from  figure  4(b)  that  the 
finite  strain  model  has  a  wave  speed  slower  than  that  in  the 
full  nonlinear  model,  and  that  waves  in  both  nonlinear  models 
are  slower  than  that  in  the  corresponding  linear  strain  model. 
Whether  this  predicted  effect  on  the  wave  speed  occurs  at  finite 
strains  remains  to  be  experimentally  validated. 

8.  Conclusions 

In  this  paper,  we  have  derived  the  three-dimensional  constitu¬ 
tive  equations  for  second-order  nonlinear  piezoelectricity  with 
6 mm  crystal  class  and  m3  crystal  class.  As  an  application, 
the  mathematical  models  of  one-dimensional  impact  problems 
are  also  completely  described.  The  one-dimensional  nonlinear 
equations  were  then  solved  numerically  via  COMSOL  and 
a  significant  effect  of  the  piezoelectric  nonlinearity  on  the 
phase  and  amplitude  of  the  stress  wave  was  clearly  observed. 
Therefore,  this  work  forms  the  basis  for  the  development  of 
nonlinear  finite  element  codes  for  the  analysis  of  dynamically 
loaded  structures  composed  of  piezoelectric  ceramic  with  the 
symmetry  class  6 mm  and/or  symmetry  class  m3. 
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